1.1 Packages nécessaires

# Librairies utilisées
library(sf)
library(cartography)
library(mapsf)
# library(corrplot)
# library(cowplot)
# library(MTA)
library(readxl)
library(ggplot2)
library(FactoMineR)
library(factoextra)
library(cluster)
library(reshape)
library(reshape2)
# library(flows)
# NB : Pour le package flows, la version la plus récente est disponible ici :
# remotes::install_github("rCarto/flows") ou # install.packages("mapsf")
# Pour obtenir une version plus ancienne (celle utilisée ici) : https://cran.r-project.org/src/contrib/Archive/flows/
# install.packages("packages/flows_1.1.1.tar.gz", repos=NULL, type="source")
# library(sp)
library(knitr)
# library(condformat) # https://cran.r-project.org/web/packages/condformat/vignettes/introduction.html
# library(units)
# library(stringr)
# library(dplyr)
# library(questionr)
# library(spdep) # Pour les matrices de contiguïté
library(rgeoda) # Pour les matrices de contiguïté
# library(rgdal) # Pour lire fichiers geodatabases

# Liste pour installer les packages si besoin :
# sf cartography mapsf readxl foreign dplyr flextable knitr stringr units condformat forcats ggplot2 rstatix questionr corrplot gtsummary broom GGally effects forestmodel ggeffects labelled cowplot spdep rgeoda

1.2 Import des couches géographiques consolidées

Les données utilisées sont une compilation de données principalement socio-économiques, localisées. La source principale est l’INSEE. Pour plus de détails, cf. le Markdown “Communes nouvelles : préparation des données”. NB : Ont été exclues des analyses l’Outre-mer et la Corse, espaces non concernés par les communes nouvelles.

Les géométries sont ici importées.

geom2011 <- st_read("data/geom.gpkg", layer = "geom2011", quiet = TRUE) 
geom_new <- st_read("data/geom.gpkg", layer = "geom_new", quiet = TRUE) 
geomfus2011 <- st_read("data/geom.gpkg", layer = "geomfus2011", quiet = TRUE) 
geomCN_new <- st_read("data/geom.gpkg", layer = "geomCN_new", quiet = TRUE)  
dep <- st_read("data/geom.gpkg", layer = "dep", quiet = TRUE)

# Pour faciliter la cortographie
# ShpDep <- as(dep, "Spatial")
ShpDep <- as_Spatial(dep) # Pose moins de problème

# Métadonnées
# Liste toutes les variables disponiles
# variables_dispo <- as.data.frame(read_excel("data-raw/meta.xlsx", sheet = "ind_target"))
# Liste les variables marquées dans le fichier meta_budgets.xlsx comme nous intéressant
# target <- subset(variables_dispo, variable_selec == "X")

1.3 Import des géométries des arrondissements

# Import des géométries des arrondissements
# Données transmises par Ludovic Chalonge (suite à une demande de Nicolas Verdier), le 30 août 2023.

# ## Ancienne version en utilisant le package rgdal
# The input file geodatabase
# fgdb <- "data-raw/arrdtv3.gdb"
# # List all feature classes in a file geodatabase
# subset(ogrDrivers(), grepl("GDB", name))
# fc_list <- ogrListLayers(fgdb)
# print(fc_list)
# date_choisie <- "1931"
# # Read the feature class
# arrondissements <- readOGR(dsn=fgdb,layer= paste0("aarrdt", date_choisie))
# projCodes <- make_EPSG() # Lister les différentes projections existant
# selecRgf <- projCodes[grep("RGF", projCodes$note), ] # Sélectionner certaines projections
# ShpDep2 <- spTransform(x = ShpDep, CRSobj = CRS("+init=epsg:2154")) # Transformer un objet dans la projection RGF

st_layers(dsn = "data-raw/arrdtv3.gdb")
date_choisie <- "1931"
arrondissements <- st_read("data-raw/arrdtv3.gdb", layer= paste0("aarrdt", date_choisie))
rm (date_choisie)


plot(arrondissements[1])

arrondissements <- st_transform(x = arrondissements, crs = st_crs(geom2011))

Les données socio-économiques qui décrivent les communes en géographies 2011 et 2021 sont ici importées. On commence par extraire les communes ayant participé à la création d’une commune nouvelle, appelées ici communes fusionnantes (datafus2011), les communes nouvelles, avec les géométries au 1er janvier 2021 et caractérisées par les données à la géométrie 2011 agrégées (dataCN_new), ainsi que les communes, à la géométrie 2011, qui n’ont pas participé à la création d’une commune nouvelle (dataNfus2011)

load("data/refdata.Rdata")
# datafus2011 <- subset(df2011, COM_NOUV == "OUI")
# dataCN_new <- subset(df_new, COM_NOUV == "OUI")
# dataNfus2011 <- subset(df2011, COM_NOUV == "NON") 

1.4 Jointure données / géométries

Dans un certain nombre de cas, il sera utile d’avoir, dans un même objet, les données et les géométries. Les données sont ici jointes aux couches géographiques d’intérêt.

geom2011 <- merge(geom2011, df2011, by = "CODGEO")
geom_new <- merge(geom_new, df_new, by = "CODGEO_new")
# geomCN_new <- merge(geomCN_new, dataCN_new, by = "CODGEO_new")
# geomfus2011 <- merge(geomfus2011, datafus2011, by = "CODGEO")

1.5 Définition de sous-ensembles pour les tests

Un certain nombre d’espaces peuvent être étudiés à part : certaines communes nouvelles car elles ont été l’objet d’une enquête qualitative (Bideau (2017), Bideau (2019)), d’autres car elles sont sur des territoires très touchés par les fusions communales. Ces sous-ensembles sont cartographiés dans la section 1.7.

testEdC <- subset (geom2011, CODGEO_new == "61324" | CODGEO_new =="73150" | CODGEO_new == "73006")
testEdCSavoie <- subset (geom2011, CODGEO_new =="73150" | CODGEO_new == "73006")
testNormandie <- subset (geom2011, REG == "23" | REG == "25") # Haute et Basse Normandie
test49 <- subset (geom2011, CODE_DEPT == "49" ) # Département du Maine-et-Loire
testOuest <- subset (geom2011, REG == "23" | REG == "25"| REG == "53"| REG == "52") # Normandies, Bretagne, Pays-de-la-Loire

1.6 Une base très riche : la base Cassini

Le site Cassini a été créé par l’EHESS, permet de retrouver des évolutions très nombreuses concernant les communes françaises.

Lien donné par Nicolas Verdier (juillet 2023) : https://didomena.ehess.fr/concern/data_sets/6395wb092?locale=fr (base en cours d’attribution d’un DOI, possiblement réalisation d’un data paper).

# download.file("https://didomena.ehess.fr/downloads/kd17cw827?locale=fr", "data-raw/temps_long/lieux_cassini_devenus_communes.csv")
# Téléchargement effectué le 6 mars 2023

cassini <- read.csv("data-raw/temps_long/lieux_cassini_devenus_communes.csv",
                            sep=",", colClasses = "character", head = TRUE, stringsAsFactors = TRUE, fill = TRUE)


# Pour simplifier la suite, on remplace dans la base Cassini les mentions "an3" et "an8" par "1795" et "1800".
# Syntaxe qui nécessite de vérifier que c'est la bonne colonne qu'on modifie :
# colnames(cassini)[30]
colnames(cassini)[30] <- "district_1795"
# colnames(cassini)[31]
colnames(cassini)[31] <- "departement_1795"
# colnames(cassini)[35]
colnames(cassini)[35] <- "pop_1795_info"
# colnames(cassini)[36]
colnames(cassini)[36] <- "pop_1795_val"
# colnames(cassini)[37]
colnames(cassini)[37] <- "pop_1800_info"
# colnames(cassini)[38]
colnames(cassini)[38] <- "pop_1800_val"

# Tentative de liaison entre la base cassini et la base DAC
load("Archives/data_prep 2011-2021(01)/data/refdata.Rdata")

cassini <- merge(cassini, df_new[, c("CODGEO_new", "NbrComFus")], by.x = "commune_mars_2021", by.y = "CODGEO_new")
# Attention, dans la base Cassini, les communes_mars_2021 sont présentes plusieurs fois

# La variable "departement_1999" permet de localiser l'ensemble des communes. Ce qui permet d'isoler certaines parties du territoire
# table(is.na(cassini$commune_mars_2021))

# Pour la commune de Vornay (Cher), la base indique "s" comme département. On remplace par "18".
# cassini$nom_1999[which(cassini$departement_1999 == "s")]
cassini$departement_1999[which(cassini$departement_1999 == "s")] <- 18
# On corrige simplement les départements avec un seul chiffre (1, 2, 3 etc) pour homogénéiser (transformés en 01, 02, 03 etc)
cassini$departement_1999 <- stringr::str_pad(cassini$departement_1999, 2, pad = "0")
print ("Les départements 1999 représentés dans la base Cassini")

## [1] "Les départements 1999 représentés dans la base Cassini"

table(cassini$departement_1999)

## 
##   01   02   03   04   05   06   07   08   09   10   11   12   13   14   15   16 
##  510  911  441  269  201  166  386  559  477  504  506  678  126  929  290  459 
##   17   18   19   21   22   23   24   25   26   27   28   29   30   31   32   33 
##  544  330  339  759  401  408  721  703  419  872  473  310  398  671  697  613 
##   34   35   36   37   38   39   40   41   42   43   44   45   46   47   48   49 
##  389  375  297  325  618  765  393  332  368  304  231  372  472  448  206  428 
##   50   51   52   53   54   55   56   57   58   59   60   61   62   63   64   65 
##  695  723  559  294  679  605  268  994  403  705  769  634 1022  502  683  517 
##   66   67   68   69   70   71   72   73   74   75   76   77   78   79   80   81 
##  279  593  400  318  678  643  421  344  337   15 1001  593  278  373 1017  397 
##   82   83   84   85   86   87   88   89   90   91   92   93   94   95 
##  272  163  168  336  356  253  569  498  108  210   36   40   49  199

# Dans certains enregistrements, il y a la mention "lac" alors que normalement c'est "lac.". On procède donc au remplacement
table(cassini$pop_1820_info) # Par exemple

## 
##  adm.  ill.   lac  lac.   pop    si 
##   388    37     8    32 38653  4271

cassini[cassini == "lac"] <- "lac."

# Possible de retirer certains départements
# cassini_reduite <- subset(cassini, departement_1999 != "2A" & departement_1999 != "2B" & departement_1999 != "73" & departement_1999 != "74")
# table(cassini_reduite$departement_1999)

1.7 Jointure base Cassini avec géométries communales

# Tranformation du data.frame cassini en objet spatial grâce au champs "position_cassini"
cassini_pos <- st_sf(
  cassini,
  geometry = st_as_sfc(cassini$position_cassini),
  crs = "EPSG:2154"
)

# Changement de la projection pour correspondre à celle des communes
cassini_pos <- st_transform(x = cassini_pos, crs = st_crs(geom2011))

# Jointure spatiale 
# Nous voulons obtenir les informations du polygones dans lequel se trouve chaque point de la géométrie cassini_pos
cassini_pos <- st_join(x = cassini_pos, y = geom2011[, "CODGEO"], join = st_intersects)
# Les dernières colonnes de cassini sont les variables (1 seule ici, CODGEO) de la table des polygones
head(cassini_pos[, c(1:4, ncol(cassini_pos)-1)])

## Simple feature collection with 6 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 3926010 ymin: 2549664 xmax: 3966608 ymax: 2573801
## Projected CRS: ETRS89-extended / LAEA Europe
##   commune_mars_2021   gid noacass featuretype CODGEO                geometry
## 1             01001 25046   60010     clocher  01001 POINT (3928662 2570876)
## 2             01001 25045      40     clocher  01001 POINT (3928602 2573801)
## 3             01002 25047   18174     clocher  01002 POINT (3966608 2554822)
## 4             01004 25048     592       bourg  01004 POINT (3961238 2549664)
## 5             01004 25049   60001      hameau  01004 POINT (3964688 2550963)
## 6             01005 25050     591     clocher  01005 POINT (3926010 2555996)

# On renome la colonne CODGEO pour davantage de clarté
colnames(cassini_pos)[ncol(cassini_pos)-1]

## [1] "CODGEO"

colnames(cassini_pos)[ncol(cassini_pos)-1] <- "CODGEO_2011"


# Problème : il y a plusieurs communes dont on ne repère pas la position ainsi
cassini_prbl_pos <- subset (cassini_pos, is.na(cassini_pos$CODGEO))
# On observe par départements
table(cassini_prbl_pos$departement_1999)

## 
## 08 13 17 22 25 29 33 34 44 50 56 57 59 64 66 68 74 80 83 85 
##  5  6  9  2  2  6  1  1  3  1  3 14  4  1  2  1  3  2  6  2

# Et cartographiquement
plot(cassini_prbl_pos$geometry)

nrow(cassini_prbl_pos)

## [1] 74

# La très grande majorité de ces communes sont sur les frontières françaises, mais elles ne sont que 74.
# Dans le cas où une part sont liées à des questions sur la Corse
# nrow (subset(cassini_prbl_pos, cassini_prbl_pos$departement_1999 != "2A" & cassini_prbl_pos$departement_1999 != "2B"))


rm(cassini_prbl_pos)

1.8 Ajout des données de population 2011 et récentes [Jachère]

# On joint les données de population à la géométrie 2011
cassini_test <- merge (cassini_pos, df2011[, c("CODGEO", "P09_POP")], by.x = "CODGEO_2011", by.y = "CODGEO")
# On renome pour plus de clarté
colnames(cassini_test)[ncol(cassini_test)-1]
colnames(cassini_test)[ncol(cassini_test)-1] <- "pop_2011_val"
# N'EST PAS ADAPTÉ PUISQUE LES POPULATIONS SONT PRÉSENTES POUR TOUTES LES COMMUNES ET PAS SEULEMENT POUR CELLES QUI SONT ENCORE ACTIVES

1.9 Jointure Cassini - Arrondissements [Jachère]

Laissé de côté pour l’instant.

# Jointure spatiale 
# Nous voulons obtenir les informations du polygones dans lequel se trouve chaque point de la géométrie cassini_pos
cassini_pos <- st_join(x = cassini_pos, y = arrondissements[, 1], join = st_intersects)

# Les dernières colonnes de cassini sont les variables (1 seule ici, CODGEO) de la table des polygones
head(cassini_pos[, c(1:4, ncol(cassini_pos)-1)])

# On renome la colonne CODGEO pour davantage de clarté
# colnames(cassini_pos)[ncol(cassini_pos)-1]
# colnames(cassini_pos)[ncol(cassini_pos)-1] <- "nom_arrond"

2.1 Évolution sommaire

# 1793-1999
# Source : http://cassini.ehess.fr/fr/html/8_glossaire.php?d=BDCarto à la définition du nom "commune"
annees <- c(1793, 1800, 1806, 1876, 1936, 1968, 1975, 1999)
nbr_com <- c(40431, 40345, 40210, 37745, 38014, 37708, 36296, 36565)
# 1876  y compris l’Alsace-Lorraine
graph <- data.frame(annees, nbr_com)

# ajout données 2023 pour voir (sources INSEE, France métropolitaine + oute-mer)
graph[nrow(graph)+1, ] <- c(1982, 36433)
graph[nrow(graph)+1, ] <- c(2012,   36571)
graph[nrow(graph)+1, ] <- c(2023, 34945)


# plot(annees, nbr_com)

ggplot(data=graph, aes(x=annees, y=nbr_com)) +
  geom_line() +
  labs(title="Évolution du nombre de communes en France", subtitle="1793-2023", x="", y = "", caption="Source : cassini.ehess.fr")

ggplot(data=graph, aes(x=annees, y=nbr_com)) +
  geom_line() +
  ylim(0, 41000) +
  labs(title="Évolution du nombre de communes en France", subtitle="1793-2023", x="", y = "", caption="Source : cassini.ehess.fr")

rm (annees)

2.2 Exploration graphique de la base Cassini

Proposition d’une exploitation pour avoir quelques graphiques un peu propres.

# On peut aussi utiliser la colonne "intendances_types"
print ("Les types d'intendances représentés dans la base Cassini")

## [1] "Les types d'intendances représentés dans la base Cassini"

table(cassini$intendances_types)

## 
##                                               
##                                             2 
##                                    intendance 
##                                         41322 
##                         intendance,intendance 
##                                           962 
##              intendance,intendance,intendance 
##                                            21 
##             intendance,souveraineté étrangère 
##                                            16 
##                        souveraineté étrangère 
##                                          1025 
##             souveraineté étrangère,intendance 
##                                            26 
## souveraineté partagée ou contestée,intendance 
##                                            15

# On décide d'enlever les territoires qui, au départ, étaient au moins en partie sous souveraineté étrangère
cassini_reduite <- subset(cassini, intendances_types != "intendance,souveraineté étrangère" & intendances_types != "souveraineté étrangère" & intendances_types != "souveraineté étrangère,intendance" & intendances_types != "souveraineté partagée ou contestée,intendance")
# On enlève également les territoires perdus en 1870 (Alsace, Moselle)
cassini_reduite <- subset(cassini_reduite, departement_1999 != "57" & departement_1999 != "67" & departement_1999 != "68")
# table(cassini_reduite$departement_1999)

# Test pour ne sélectionner que l'Alsace-Moselle
# cassini_reduite <- subset(cassini, departement_1999 == "57" | departement_1999 == "67" | departement_1999 == "68")



# En l'absence de méta-données détaillées, on choisit de ne s'intéresser qu'aux communes dont la population est renseignée ou simplement lacunaire (devrait être renseignée). On postule que les autres indications désignent des communes n'existant pas à l'époque.
print ("Les types d'informations renseignées pour le recensement de l'an 3 (1794-1795)")

## [1] "Les types d'informations renseignées pour le recensement de l'an 3 (1794-1795)"

table (cassini_reduite$pop_1795_info)

## 
##  abs.  adm.  lac.   pop    si 
##   215    31   136 36917  3102

# On cherche les noms des colonnes qui contiennent le motif "_info"
labels_data_pop <- colnames(cassini[grep(pattern = "_info", colnames(cassini), ignore.case = FALSE)])
label <- labels_data_pop [2] # Si besoin pour tests

evol_nbr_com <- data.frame() # Création d'un objet pour récupérer les résultats totaux
evol_nbr_com_reduit <- data.frame() # Création d'un objet pour récupérer les résultats de la base réduite

# Boucle pour exploiter les données date par date
for (label in labels_data_pop) {
  annee <- substr(label, start =5, stop = nchar(label)-5) # On repère le n° de l'année en segmentant le label
  sous_tableau <- subset(cassini, cassini[, label] == "pop" | cassini[, label] == "lac." | cassini_reduite[, label] == "abs.") # On crée un sous-tableau comptant uniquement les communes dont la population nous intéresse, à l'année donnée
  # sous_tableau <- subset(cassini, cassini[, label] == "pop") # Juste avec les populations connues
  nbr_com <- nrow(sous_tableau) # On compte le nombre de communes
  
  sous_tableau_reduit <- subset(cassini_reduite, cassini_reduite[, label] == "pop" | cassini_reduite[, label] == "lac." | cassini_reduite[, label] == "abs.") # On crée un sous-tableau comptant uniquement les communes dont la population nous intéresse, à l'année donnée
  # sous_tableau_reduit <- subset(cassini_reduite, cassini_reduite[, label] == "pop") # Juste avec les populations connues
  nbr_com_reduit <- nrow(sous_tableau_reduit) # On compte le nombre de communes
  
  i <- c(annee, nbr_com) # On compile les résultats totaux
  i_reduit <- c(annee, nbr_com_reduit) # On compile les résultats totaux
  evol_nbr_com <- rbind(evol_nbr_com, i, stringsAsFactors = FALSE)
  evol_nbr_com_reduit <- rbind(evol_nbr_com_reduit, i_reduit, stringsAsFactors = FALSE)
}

# On traite à la main les deux années 2011 et 2021
# Pour 2011, on fait la liste des CODGEO unique de 2021 et on multiplie par leur nombre de communes fusionnantes (donne le nombre de communes en 2011) 
# Pas très propre mais pas l'énergie de faire mieux là...
CODGEO_2021 <- as.data.frame(unique(cassini$commune_mars_2021))
colnames(CODGEO_2021) <- "CODGEO_new"
CODGEO_2021_reduit <- as.data.frame(unique(cassini_reduite$commune_mars_2021))
colnames(CODGEO_2021_reduit) <- "CODGEO_new"
# On joint pour avoir le nombre de communes fusionnantes
CODGEO_2021 <- merge(df_new[c("CODGEO_new", "NbrComFus")], CODGEO_2021, by = "CODGEO_new", all.y = TRUE)
CODGEO_2021_reduit <- merge(df_new[c("CODGEO_new", "NbrComFus")], CODGEO_2021_reduit, by = "CODGEO_new", all.y = TRUE)
# Pour les communes non fusionnantes, on indique qu'il y en a une
CODGEO_2021$NbrComFus[is.na(CODGEO_2021$NbrComFus)] <- 1
CODGEO_2021_reduit$NbrComFus[is.na(CODGEO_2021_reduit$NbrComFus)] <- 1

i_2011 <- c(2011, sum(CODGEO_2021$NbrComFus))
i_2011_reduit <- c(2011, sum(CODGEO_2021_reduit$NbrComFus))
evol_nbr_com <- rbind(evol_nbr_com, i_2011)
evol_nbr_com_reduit <- rbind(evol_nbr_com_reduit, i_2011_reduit)


# Pour 2021, on liste simplement le nombre de CODGEO uniques de 2021 pour chacun des groupe
i_2021 <- c(2021, length(unique(cassini$commune_mars_2021)))
i_2021_reduit <- c(2021, length(unique(cassini_reduite$commune_mars_2021)))

evol_nbr_com <- rbind(evol_nbr_com, i_2021)
evol_nbr_com_reduit <- rbind(evol_nbr_com_reduit, i_2021_reduit)

colnames (evol_nbr_com) <- c("annees", "nbr_com")
colnames (evol_nbr_com_reduit) <- c("annees", "nbr_com")
evol_nbr_com$base <- "tout"
evol_nbr_com_reduit$base <- "tjr_fr"

evol_nbr_com <- rbind(evol_nbr_com, evol_nbr_com_reduit)

# evol_nbr_com[c(1,2), 1] <- c(1795, 1800) # On modifie l'intitulé pour faciliter la représentation graphique
evol_nbr_com$annees[evol_nbr_com$annees == "an3"] <- 1795
evol_nbr_com$annees[evol_nbr_com$annees == "an8"] <- 1800

evol_nbr_com[,1] <- as.numeric(evol_nbr_com[,1]) # On passe en numérique pour faciliter la représentation
evol_nbr_com[,2] <- as.numeric(evol_nbr_com[,2]) # Si on veut faire une seule colonne



ggplot(data=evol_nbr_com, aes(x=annees, y=nbr_com, group=base)) +
  geom_line(aes(group = base, linetype = base)) +
  scale_linetype(name = "Périmètre", labels = c("Territoires toujours français", "Toute base Cassini")) +
  scale_x_continuous (breaks = seq(1780, 2030, 20)) +
  labs(title="Évolution du nombre de communes", subtitle="1795-2021", x="", y = "", caption="Sources : cassini.ehess.fr, INSEE")

ggplot(data=evol_nbr_com, aes(x=annees, y=nbr_com, group=base)) +
  geom_line(aes(group = base, linetype = base)) +
  scale_linetype(name = "Périmètre", labels = c("Territoires toujours français", "Toute base Cassini")) +
  ylim(0, 41000) +
  labs(title="Évolution du nombre de communes", subtitle="1795-2021", x="", y = "", caption="Sources : cassini.ehess.fr, INSEE")

rm(label, annee, sous_tableau, nbr_com, i, nbr_com_reduit, sous_tableau_reduit, i_reduit, i_2021, i_2021_reduit, CODGEO_2021, CODGEO_2021_reduit, i_2011, i_2011_reduit)

2.3 Se concentrer sur les données les plus fiables

La base Cassini est très riche, mais elle reste tributaire des sources, parfois peu fiables. Ainsi, les données indiquées comme “An 3” (1794-1795) peuvent être considérées comme peu fiables. Nous les retirons donc de l’analyse.

cassini <- subset(cassini, select= -c(pop_1795_info, pop_1795_val))
labels_data_pop <- colnames(cassini[grep(pattern = "_info", colnames(cassini), ignore.case = FALSE)])

4.1 Réalisation d’un graphique par département

liste_dep <- sort(unique (cassini$departement_1999))

num_dep <- 21

# Création d'un tableau pour accueillir les données de l'ensemble des départements
evol_deps <- data.frame()

for (num_dep in liste_dep) {
  
cassini_reduite <- subset(cassini, departement_1999 == num_dep)
nom_dep <- dep$LIBELLE[dep$CODE_DEPT == num_dep]


evol_nbr_com_reduit <- data.frame() # Création d'un objet pour récupérer les résultats de la base réduite

# label <- labels_data_pop[2]

# Boucle pour exploiter les données date par date
for (label in labels_data_pop) {
  annee <- substr(label, start =5, stop = nchar(label)-5) # On repère le n° de l'année en segmentant le label
  sous_tableau_reduit <- subset(cassini_reduite, cassini_reduite[, label] == "pop" | cassini_reduite[, label] == "lac." | cassini_reduite[, label] == "abs.") # On crée un sous-tableau comptant uniquement les communes dont la population nous intéresse, à l'année donnée
  # sous_tableau_reduit <- subset(cassini_reduite, cassini_reduite[, label] == "pop") # Juste avec les populations connues
  nbr_com_reduit <- nrow(sous_tableau_reduit) # On compte le nombre de communes
  pop_annee_dep <- sum(as.numeric(sous_tableau_reduit[, paste0("pop_", annee, "_val")]), na.rm=TRUE)
  i_reduit <- c(annee, nbr_com_reduit, pop_annee_dep) # On compile les résultats totaux
  evol_nbr_com_reduit <- rbind(evol_nbr_com_reduit, i_reduit, stringsAsFactors = FALSE)
}

# On traite à la main les deux années 2011 et 2021
# Pour 2011, on fait la liste des CODGEO unique de 2021 et on multiplie par leur nombre de communes fusionnantes (donne le nombre de communes en 2011) 
# Pas très propre mais pas l'énergie de faire mieux là...
CODGEO_2021_reduit <- as.data.frame(unique(cassini_reduite$commune_mars_2021))
colnames(CODGEO_2021_reduit) <- "CODGEO_new"
# On joint pour avoir le nombre de communes fusionnantes
CODGEO_2021_reduit <- merge(df_new[c("CODGEO_new", "NbrComFus", "P09_POP")], CODGEO_2021_reduit, by = "CODGEO_new", all.y = TRUE)
# Pour les communes non fusionnantes, on indique qu'il y en a une
CODGEO_2021_reduit$NbrComFus[is.na(CODGEO_2021_reduit$NbrComFus)] <- 1

i_2011_reduit <- c(2011, sum(CODGEO_2021_reduit$NbrComFus), sum(CODGEO_2021_reduit$P09_POP))
evol_nbr_com_reduit <- rbind(evol_nbr_com_reduit, i_2011_reduit)


# Pour 2021, on liste simplement le nombre de CODGEO uniques de 2021 pour chacun des groupe
i_2021_reduit <- c(2021, length(unique(cassini_reduite$commune_mars_2021)), NA)

evol_nbr_com_reduit <- rbind(evol_nbr_com_reduit, i_2021_reduit)

colnames (evol_nbr_com_reduit) <- c("annees", "nbr_com", "population")

evol_nbr_com <- evol_nbr_com_reduit

# evol_nbr_com[c(1,2), 1] <- c(1795, 1800) # On modifie l'intitulé pour faciliter la représentation graphique
# Est désormais fait plus haut
# evol_nbr_com$annees[evol_nbr_com$annees == "an3"] <- 1795 
# evol_nbr_com$annees[evol_nbr_com$annees == "an8"] <- 1800

evol_nbr_com[,1] <- as.numeric(evol_nbr_com[,1]) # On passe en numérique pour faciliter la représentation
evol_nbr_com[,2] <- as.numeric(evol_nbr_com[,2]) # Si on veut faire une seule colonne
evol_nbr_com[,3] <- as.numeric(evol_nbr_com[,3]) # Si on veut faire une seule colonne


# Calcul de l'évolution en valeur absolue du nombre de commune dans le département
evol_nbr_com$evol <- evol_nbr_com$nbr_com - lag(evol_nbr_com$nbr_com)

# Calcul du taux d'évolution
evol_nbr_com$tx_evol <- round(100 * (evol_nbr_com$nbr_com - lag(evol_nbr_com$nbr_com)) / lag(evol_nbr_com$nbr_com), 2)

# On compile les résultats du département en question
departement <- c(nom_dep, num_dep, evol_nbr_com$nbr_com, evol_nbr_com$population, evol_nbr_com$evol, evol_nbr_com$tx_evol)
# On ajoute ces résultats dans le tableau global
evol_deps <- rbind (evol_deps, departement, stringsAsFactors = FALSE)

# ggplot(data=evol_nbr_com, aes(x=annees, y=nbr_com)) +
#   geom_line() +
#   scale_x_continuous (breaks = seq(1780, 2030, 20)) +
#   labs(title="Évolution du nombre de communes en Aveyron", subtitle="1795-2021", x="", y = "", caption="Sources : cassini.ehess.fr, INSEE")

graphique <- ggplot(data=evol_nbr_com, aes(x=annees, y=nbr_com)) +
  geom_line() +
  ylim(0, max(evol_nbr_com$nbr_com)+50) +
  scale_x_continuous (breaks = seq(1780, 2030, 20)) +
  labs(title= paste0("Évolution du nombre de communes en ", nom_dep, " (", num_dep, ")"), subtitle="1795-2021", x="", y = "", caption="Sources : cassini.ehess.fr, INSEE")

print(graphique)

# Pour avoir un graphique plus complet, on change la forme des données
pr_graph <- melt(evol_nbr_com, id.vars = "annees")
  

ggplot(data=pr_graph, aes(x=annees, y = value)) +
  geom_line(aes(color = variable, linetype = variable)) +
  # ylim(min(evol_nbr_com$evol), max(evol_nbr_com$nbr_com)+50) +
  scale_x_continuous (breaks = seq(1780, 2030, 20)) + 
  scale_linetype(name = "Évolution", labels = c("Nombre de communes", "Différence entre recensements", "Taux d'évolution entre recensements")) +
  scale_color_discrete(name = "Évolution", labels = c("Nombre de communes", "Différence entre recensements", "Taux d'évolution entre recensements")) +
  labs(title= paste0("Communes par départements : ", nom_dep, " (", num_dep, ")"), subtitle="1795-2021", x="", y = "", caption="Sources : cassini.ehess.fr, INSEE")

rm(label, annee, nbr_com_reduit, sous_tableau_reduit, i_reduit, i_2021_reduit, CODGEO_2021_reduit, i_2011_reduit, graphique)

}

# Création d'une colonne identifiant les deux recensements en jeu dans les évolutions
evol_nbr_com$ident_evol <- paste0("evol_", lag(evol_nbr_com$annees), "_", evol_nbr_com$annees)

# On identifie les colonnes du tableau regroupant les données de tous les départements
colnames(evol_deps) <- c("Département", "CODE_DEPT", evol_nbr_com$annees, paste0("pop_", evol_nbr_com$annees), evol_nbr_com$ident_evol, paste0("tx_", evol_nbr_com$ident_evol))

# print (evol_deps)
evol_deps_affichage <- evol_deps
colnames(evol_deps_affichage) <- c("Département", "N° Département", evol_nbr_com$annees, paste0("pop_", evol_nbr_com$annees), evol_nbr_com$ident_evol, paste0("tx_", evol_nbr_com$ident_evol))

nbr_col <- length(labels_data_pop) + 2

kable(evol_deps_affichage, row.names = F, digits = 2, align = "c",
      caption = "Évolution des communes par départements (à la géographie 1999)") %>%
  kableExtra::add_header_above(c(" " = 2, "Nombre de communes" = nbr_col, "Population" = nbr_col, "Évolution en valeur absolue" = nbr_col, "Taux d'évolution entre chaque recensement" = nbr_col))